Defect Modes in One-Dimensional Granular Crystals 
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We study the vibrational spectra of one-dimensional statically compressed granular crystals (ar- 
rays of elastic particles in contact) containing defects. We focus on the prototypical settings of one 
or two spherical defects (particles of smaller radii) interspersed in a chain of larger uniform spherical 
particles. We measure the near-linear frequency spectrum within the spatial vicinity of the defects, 
and identify the frequencies of the localized defect modes. We compare the experimentally deter- 
mined frequencies with those obtained by numerical eigen-analysis and by analytical expressions 
based on few-site considerations. We also present a brief numerical and experimental example of 
the nonlinear generalization of a single-defect localized mode. 

PACS numbers: 63.20.Ry, 63.20.Pw, 45.70.-n, 46.40.-f 



INTRODUCTION 

Defect modes in crystals have long been studied in the 
realm of solid state physics P, Q • The presence of defects 
or "disorder" is known to enable localized lattice vibra- 
tions, whose associated frequencies have been measured 
in the spectra of real crystals (see 0, Q and references 
therein). More recently, this study has been extended 
to include other examples, including superconductors Q 
and electron-phonon interactions [5] . Similar phenomena 
have also been observed in nonlinear systems, includ- 
ing photonic crystals [H, optical waveguide arrays @- 
9[_, dielectric superlattices (with embedded defect layers) 
[lfj , micromechanical cantilever arrays and Bose- 
Einstein condensates of atomic vapors [12, LL3J ■ 

Granular crystals are nonlinear systems composed of 
densely-packed particles interacting through Hertzian 
contacts [l4T - [l7j ]. These systems present a remarkable 
ability to tune their dynamic response from linear to 
strongly nonlinear regimes [H|. This has allowed the 
exploration of fundamental nonlinear waveforms such as 
traveling waves 15 - 3] and discrete breathers [l9| . Gran- 
ular crystals have also been proposed for several e ngineer 
ing applications, such as energy absorbing layers [2J 
actuating devices 24 1, and sound scramblers [25II2 



The presence of defects in statically uncompressed (or 
weakly compressed, as compared to the relative dynamic 
displacements) granular chains excited by impulsive load- 
ing has been studied in a number of previous works that 
have reported the existence of interesting dynamic re- 
sponses such as the fragmentation of waves, anomalous 



reflections, and energy trapping 2pl- 23l l27H32l |. In this 
paper, we study the response of strongly compressed 
granular crystals, with one or two defects (extending our 
earlier theoretical work 33[), excited by continuous sig- 
nals. We measure the frequency response of the system 
and reveal localized modes due to the presence of de- 
fects. We report that the number of localized modes 
mirrors that of the defects, and note that the frequen- 



cies of such modes depend on: (i) the ratio of the defect 
mass to the mass of the particles in the uniform chain, 
(ii) the relative proximity of multiple defects, (iii) the 
geometric and material properties of the particles com- 
posing the crystal, and (iv) the static load. We compare 
our experimental findings with numerical computations 
and with theoretical analysis approximating the behavior 
of a few sites in the vicinity of the defect (s). Finally, we 
demonstrate that as we go from the linear to the nonlin- 
ear regime, nonlinear "deformations" of the linear defect 
modes (with appropriately downshifted frequencies) are 
sustained by the system. 



EXPERIMENTAL SETUP 

We assemble I D g ranular crystals, similar to those 
described in [l9|, |34| - composed of N — 20 statically 
compressed stainless steel spherical particles (316 type, 
with elastic modulus E = 193 GPa and Poisson ra- 
tio Vh — 0.3 [HI), as shown in Fig. QJa). The chain 
is composed of uniform particles of (measured) radius 
R = 9.53 mm and mass M = 28.84 g, except for one 
(or two) light-mass stainless steel defect particles. The 
spheres are held in a ID configuration using four polycar- 
bonate bars (12.7 mm diameter) that are aligned by poly- 
carbonate guide plates spaced at approximately 12 cm in- 
tervals along the axis of the crystal. The defect particles, 
which are of smaller radii than the rest of the particles of 
the chain, are aligned with the axis of the crystal using 
polycarbonate support rings. Dynamic perturbations are 
applied to the chain by a piezoelectric actuator mounted 
on a steel cube (which acts as a rigid wall). The parti- 
cles are statically compressed by a load of F = 20 N. 
The static load is applied using a soft spring (of stiffness 
1.24 kN/m), which is compressed between the last par- 
ticle in the chain and a second steel cube bolted to the 
table. The applied static load is measured by a calibrated 
load cell placed between the spring and the steel cube. 
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We measure the dynamic force signals of the propagat- 
ing waves with custom-made force sensors consisting of a 
piezo-electric disk embedded inside two halves of a stain- 
less steel particle with radius R = 9.53 mm. The sensor- 
particles are carefully constructed to resemble the mass, 
shape, and contact properties of the other spherical par- 
ticles composing the rest of the crystal 2(J 25, 26, 37 ,l38j| . 



displacement obtained from a static load Fq is A„ = 
(F /A n ) 2 / 3 [lj, [l5|, where R n is the radius of the nth 
particle. 

In order to study the linear spectrum of the inhomo- 
geneous granular crystal, we linearize Eq. ((T|) about the 
equilibrium state under the presence of the static load. 
This yields the following linear system [l9l . [33l [3(| : 



a). 
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m n u n = K n (u n -i - u n ) - K n+ i(u n - u n+ i) 
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FIG. 1: [Color online] a) Schematic diagram of the experi- 
mental setup for the homogeneous chain with a single defect 
configuration, b) Experimental transfer functions (as defined 
in the "single-defect: near linear regime" section) for a gran- 
ular crystal with a static load of Fo = 20 N and a defect-bead 
of mass m = 5.73 g located at site ridef = 2. Blue (dark-grey) 
[red (light-grey)] curves corresponds to transfer function ob- 
tained from the force signal of a sensor particle placed at n — 4 
[n = 20]. The diamond marker is the defect mode. The tri- 
angle marker is the upper acoustic cutoff mode. The vertical 
black dashed line is the theoretically predicted defect mode 
frequency, and the vertical solid black line is the theoretically 
predicted upper acoustic cutoff frequency. 



THEORETICAL MODEL 

We consider the ID inhomogeneous crystal of N beads 
as a chain of nonlinear oscillators 1 151: 



m n u n = A n [A n + u n -\ - u, 
- A n+ i[A n+ i +u n - 



(1) 



where [Y] + denotes the positive part of Y (which signifies 
that adjacent particles interact only when they are in 
contact), u n is the displacement of the nth sphere (where 
n£ {1, • • • , N}) around the static equilibrium, m n is the 
mass of the nth particle, and the coefficients A n depend 
on the exponent p and the geometry/material properties 
of adjacent beads. The exponent p = 3/2 represents the 
Hertz law potential between adjacent spheres [lij]. In 

1/2 



this case, A n 
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and the static 



where K n — \A 2 n Z 'F^ 3 . Following [33[, we simplify 
Eq.([2]) to the eigensystem: 



w 2 Mu = Au, 



(3) 



where M is a N x N diagonal matrix with elements 
M nn — m-m and u is the displacement vector. A is a 
N x N triagonal matrix with elements A mn — —\K n + 
(l-6 nN )K n+ i]8 mn +K n+1 (5 mn - 1 +6 mn+1 ), where 5 is the 
Kronecker delta and we consider left-fixed and right-free 
boundary conditions. The right-free boundary assump- 
tion derives from the low stiffness of the static compres- 
sion spring (Fig.rjja)) as compared to the stiffness of the 
particles in contact. 



SINGLE-DEFECT: NEAR-LINEAR REGIME 

In this section, we study ID granular crystals that are 
homogeneous except for one light-mass defect bead at 
site ridef, as shown in Fig. QJa). Solving the eigenvalue 
problem of Eq. ([3]), for such a granular crystal, we ob- 
tain the eigenfrequencies and the corresponding spatial 
profile of the modes of the system. The presence of the 
single light-mass defect generates a localized mode (see 
also 32j, |33(), centered at the defect site, which we will 
refer to as the defect mode. The defect mode amplitude 



decays exponentially away from the defect site and its 
frequency fd is such that fd 



> f c , where / c = -L^i*^ 
is the upper cutoff frequency of the acoustic band of the 

,2/3 pl/3 



homogeneous host crystal (where Krr — fA'^pQ ' ' is 
the linear stiffness of the contact between two beads with 
radius R) . The spatial profile of this mode consists of ad- 
jacent particles oscillating out of phase (see inset in Fig. 
(J2]))- As the radius of the defect bead becomes smaller, 
the difference between fd and f c becomes larger, while 
the defect mode becomes more spatially localized. We 
observe that for the granular crystals studied here, with 
radii ratios of < 0.7, the defect mode involves the 
motion of up to approximately three beads, i.e., the dis- 
placements of the beads at n > ridef + 2 and n < n^e/ — 2 
are negligible. Because in this range of radii ratio the mo- 
tion of the particles can be accurately approximated by 
three beads, we consider the particles at n = ridef ± 2 
as fixed walls, in order to find an analytical approxima- 
tion for the frequency of the defect mode. Solving for 
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the eigenfrequencies of this reduced three-bead system, 
we find that the mode corresponding to the out of phase 



motion can be analytically approximated by Eq. ((4]) 
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1 / 2K Rr M + K RR m + K Rr m + — &K Rr K RR mM + [2K Rr M + {K RR + K Rr )m] 2 
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where K Rr — %A Rr F Q is the linear stiffness of the 
contact between a defect-bead and a bead of radius R. 

We conduct experiments to identify the frequency of 
the defect mode in granular crystals with a single light- 
mass defect as shown in Fig. [Ha). We place the de- 
fect particle at site ridef = 2 (close to the actuator) 
so that the energy applied by the actuator, at the de- 
fect mode frequency, will not be completely attenuated 
by the uniform crystal, which acts as a mechanical fre- 
quency filter, before it arrives at the defect site. Be- 
cause of the localized nature of the defect mode, plac- 
ing a defect particle (of radius r < 7.14 mm) at site 
ridef = 2 or further into the chain makes nearly no dif- 
ference on the frequency of the defect mode. For in- 
stance, for a defect particle of radius r = 7.14 mm, we 
numerically calculate (using Eq. [3]) the difference in the 
defect mode frequency for the cases where a defect par- 
ticle is placed at site ridef = 2 or ridef = 10, to be 
3 Hz. Conversely, because of the presence of the fixed 
boundary and the larger localization length of the de- 
fect mode, for a defect particle of r = 8.73 mm, we cal- 
culate the difference in defect mode frequency, between 
sites ridef = 2 and ridef = 10, to be 68 Hz. The defect 
particles are stainless steel spheres of smaller radii, r = 
[3.97, 4.76, 5.56, 6.35, 7.14, 7.94] mm, and measured 
masses of m = [2.08, 3.60, 5.73, 8.54, 12.09, 16.65] g, 
respectively. We experimentally characterize the linear 
spectrum of this system by applying low amplitude (ap- 
proximately 200 mN) bandwidth limited noise (3—25 kHz 
for the two smallest defect particles, and 3—15 kHz oth- 
erwise) via the piezoelectric actuator. We calculate the 
transfer functions, specific to the sensor location, by aver- 
aging the Power Spectral Densities (PSD [39j]) of 16 force- 
time histories, measured with the embedded sensors, and 
dividing by the average PSD level in the 3 — 8 kHz range 
(corresponding to the acoustic band). We embed sensors 
in particles at sites n = 4 and n = 20. In Fig. HJb) we 
show the transfer functions for the granular crystal with 
defect radius r = 5.56 mm. The red (light-grey) and 
blue (dark-grey) curves are the transfer functions for the 
sensors at sites n — 4 and n — 20, respectively. We de- 
note the experimental cutoff frequency by the triangular 
marker (found by identifying the last peak in the acous- 
tic band) and defect frequency as the diamond marker 
on the n = 4 transfer function. The vertical lines de- 



note the theoretically determined upper cutoff frequency 
of the acoustic band and the defect frequency (Eq. [4}. 
The presence of the defect mode can be clearly identified 
in the vicinity of the defect (at n — 4) , but is not visible 
far from the defect (at n = 20). 



We repeat the process of measuring the transfer func- 
tion and identifying the defect mode frequency 16 times, 
re-assembling the crystal after each repetition. In Fig. [2] 
we plot the average frequency of the 16 experimentally 
identified defect modes as a function of the mass ratio 
(blue [dark-grey] solid line connecting the closed di- 
amonds). We also plot, for comparison, the defect fre- 
quency predicted by the analytical expression of Eq. ^ 
(green [light-grey] dashed line connecting the crosses), 
and the numerical eigenanalysis of Eq. ([3]) corresponding 
to the experimental setup (black solid line connecting the 
open diamonds). The error bars on the experimental data 
are ±2<r where a is the standard deviation of the iden- 
tified defect frequencies over the 16 repetitions. Com- 
paring the analytical three-bead approximation with the 
numerical eigenfrequencies, we find an excellent agree- 
ment for mass ratios of ■jg < 0.6. Comparing the ex- 
perimental data with the numerics, we find an upshift of 



5 — 10%, similar to the upshift observed in 19|, 34 [. For 
the r — 5.56 mm defect, the average experimental defect 
frequency is /j — 13.59 kHz and the average experi- 
mental cutoff frequency is f^ xp — 8.36 kHz. In compari- 
son, the theoretical 3-bead approximation gives a defect 
frequency of f^ bead — 12.84 kHz and the eigenproblem of 
Eq. (0 gives a defect frequency of f2' um = 12.85 kHz, 
while the analytically calculated cutoff frequency was 
f c = 8.02 kHz. 



Possible reasons for these upshifts have been identi- 
fied in 1|| 34] and the references therein, such as error 
in the material parameters, nonlinear elasticity, surface 
roughness, dissipative mechanisms and misalignment of 
the particles. We note that a systematic error in the 
measurement of the static load could also cause such an 
upshift. Nevertheless, it is clear from Fig. [2] that the func- 
tional dependence of the relevant frequencies on the mass 
ratio (of defect to regular beads) is accurately captured 
by our analytical and numerical results. 
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FIG. 2: [Color online] Frequency of the defect mode, with 
defect-bead placed at ridef — 2, as a function of mass ra- 
tio m/M. Solid blue line (dark-grey, closed diamonds) cor- 
responds to experiments, solid black line (open diamonds) 
to numerically obtained eigenfrequencies (see Eq. (|3]l). and 
green dashed line (light-grey, x markers) to the analytical pre- 
diction of the three-beads approximation (see Eq. (0). The 
error bars account for statistical errors on the measured fre- 
quencies and are ±2o\ Inset: The normalized defect mode for 
— —02 



TWO-DEFECTS: NEAR-LINEAR REGIME 



We study granular crystal configurations with two 
identical light-mass defects to better understand the ef- 
fects of increasing heterogeneity on the spectral response 
of the system. The localized mode due to the pres- 
ence of a light-mass defect (for mass ratios < 0.6) 
has a spatial localization length of about three particles 
(larger particles have a greater localization length, and 
smaller particles have a shorter localization length), as 
described in the one-defect case and shown in Fig. [5] 
We can thus expect that two light-mass defects placed 
far from each other in a granular crystal (sufficiently 
outside this localization length) would have similar fre- 
quencies and mode shapes independent of the presence 
of the other. However, as the two defect particles are 
brought closer together (within the localization length), 
each mode influences the other. For a sufficiently small 
mass ratio, this results in the creation of two defect 
modes at different frequencies; one with the defect par- 
ticles moving out of phase, and the other with the de- 
fect particles moving in phase. For the case of nearest- 
neighbor identical defects, our theoretical analysis can 
be extended by using a four-particle analogy. In this 
case, using the notation s± = K Rr (M + m) + K RR m, 
S2 = -AK Rr K RR Mm + (K RR m + K Rr (M + m)) 2 , s 3 = 
si + 2K rr M and s 4 = -A{2K rr K RR + K Rr (2K rr + 
K RR ))Mm + (2K rr M + K RR m + K Rr (M + m)) 2 , we 
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The two highest frequencies correspond to the linear 
defect mode frequencies. Naturally, this analytical ap- 
proach can be extended to more distant defects, although 
we do not present such algebraically intensive cases here. 

In Fig. |31 we show the behavior of two r = 5.56 mm 
defects in a N = 20 particle granular crystal under 
F = 20 N static load (similiar to the configuration shown 
in Fig.QJa)), where the first defect is at site k — n de fi = 2 
and the second defect is at a variable position between 
site I = ride f 2 — 3 and I = ride f 2 — 6. We use the same 
experimental method as in the single defect case except 
now we use a noise range between 3 — 20 kHz and we 
place the first sensor at n = I + 1. We show the experi- 
mentally determined PSD transfer function for the case 
of I — k = 1 in Fig. Ela), with sensors at site n = 4 (blue 
[dark-grey] ) and n = 20 (red [light-grey]). As described 



the existence of two separate defect modes for 
the case where the defect particles are adjacent to each 
other (I — k = 1), depends on the mass ratio of the de- 
fect particles to those of the rest of the crystal. Here the 
mass ratio is such that two modes are present, as can 
be seen in the blue (dark- grey) curve in Fig. [3£a). The 
two distinct modes, which we denote by the open square 
and closed circular markers, have frequencies above the 
acoustic band. The square markers denote the mode with 
defect beads moving out of phase, and the closed circu- 
lar marker corresponds to the mode with defect parti- 
cles moving in phase, as shown by the numerically calcu- 
lated eigenmodes in Fig.[3][c) and (d) respectively [HI . In 
Fig.EJb) we plot the experimentally determined frequen- 
cies of both modes as a function of the inter-defect parti- 
cle distance (I — k). The solid blue (dark-grey) lines are 
the experimental data, and the dashed black lines are the 
frequencies obtained from solving the eigenvalue problem 
of Eq. d3]) . The green (light-grey) x-markers denote the 
frequencies calculated with Eqs. ©-((I]), for the I — k = 1 
case. It is evident that the analytical results agree closely 
with the numerically calculated eigen-frequencies. The 
error bars on the experimental data correspond to the 
±2u standard deviation as calculated in the single-defect 
case. We see close qualitative agreement between the ex- 
perimental data and the numerical predictions, but also 
the same systematic upshift as observed in the single de- 
fect case and [H H3. From Fig. Efb) we can see that 
as the defects are placed three or more particles apart, 
the frequencies of the defect modes converged to approx- 
imately the same value, suggesting the defects respond 
independently of each other. 
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FIG. 3: [Color online] (a) Experimental transfer functions 
for a granular crystal with two defect-beads of mass ratio 
= 0.2 at ridef = 2 and n^ef ~ 3 (in contact). Blue (dark- 
grey) [red (light-grey)] curve corresponds to transfer function 
obtained from the force signal of a custom sensor placed at 
n — 4 [n = 20]. (b) Frequencies of the defect modes as a 
function of the distance between them. The solid line denotes 
experimental data, the dashed line the numerically obtained 
eigenfrequncies, and the x markers the frequencies from the 
analytical expresssions of Eqs. ([5])- (|6}. (c),(d) The normal- 
ized defect mode shapes corresponding to the defect modes 
identified in (a) with frequency of the same marker-type. 



SINGLE DEFECT: NONLINEAR LOCALIZED 
MODES 

As shown in 33], the interplay of the inherent non- 
linearity of the granular crystal with the linear local- 
ization due to the defect results in the presence of ro- 
bust nonlinear localized modes (NLMs). The frequency 
of these modes depends not only on the static load and 
the material values of the beads but also on the am- 
plitude of the oscillations. In order to find this depen- 
dence, we apply Netwon's method (see [33| and refer- 
ences therein) for the experimental, single-defect, con- 
figuration of Fig. 1(a). For the numerical calculations 
in this section, we calculate experimental contact coeffi- 
cients following a procedure similar to the one described 
in [l9[. The experimental contact coefficients obtained 
are A e ^ = 10.79 N//j,to 3 / 2 for the contact between two 



R = 9.53 mm beads and A c r x R p = 9.95 N/ ^m 3 / 2 for the 
contact between the R = 9.53 mm and the r — 5.56 mm 
beads. In comparison, the values of the coefficient A, as 
calculated by the material values and used for the previ- 
ous sections of the paper, are Arr = 9.76 N//zm 3 / 2 and 
A Rr = 8.38 N/ ^m 3 / 2 . 

In Fig. 4(a), we show the frequency of the numerically 
determined NLM as a function of the averaged dynamic 
force for the particle at site n = 3. The latter corresponds 
to the average of the two dynamic contact forces adja- 



cent to the particle, which is analogous to what is mea- 
sured experimentally by the dynamic force sensor [25| . 
In Fig. 4(b), we plot the numerically determined nor- 
malized NLM shape at — 13.28 kHz. Comparing this 
NLM shape to the linear mode shape of the same fre- 
quency (see inset of Fig. [5]), we can see that the NLM 
has a slightly modified (more asymmetric) spatial pro- 
file. 

The experimental setup used for the study of the NLMs 
is the same as in the case of the linear single defect exper- 
iments (as shown in Fig. [Ha)) except we place sensors in 
particles at sites n = 3, n = 5, and n = 20. Additionally, 
we replace the n — 1 particle with an embedded actuator 
particle, so as to apply high amplitude (approximately 
10 N), short time pulse (approximately 100 /xs) perturba- 
tions directly to the defect particle. Exciting such a pulse 
creates an initial condition in the crystal that resembles 
the predicted defect NLM shape. The embedded actua- 
tor particle is similar in construction to the sensors but 
with a piezoceramic construction/geometry more appro- 
priate for high force amplitude actuation (Piezomechanik 
PCh 150/5x5/2 Piezo-chip). 

The force-time history of the dynamic force measured 
by the sensor at site n = 3 is shown in Fig.^Jc). A sharp 
excitation is evident at time t = 0, followed by periodic 
oscillations with a decaying envelope, due to the inherent 
dissipation in the system. As shown by the parametric 
continuation in Fig.|U[a), NLMs corresponding to the de- 
fect mode at higher amplitudes have a frequency deeper 
into the gap than its linear counterpart. However, for the 
amplitudes observed here this is only a slight shift (up to 
200 Hz over 7 N). 

We study, in more detail, the response of two selected 
time regions of the force-time history shown in Fig. EJc) , 
to experimentally demonstrate the frequency shift char- 
acteristic of higher amplitude NLMs. The two non- 
overlapping time regions are of length T — 5.1 ms. The 
red (light-grey) time region begins immediately following 
the arrival of the initial actuated pulse, and presents a 
maximum amplitude of 7 N. The blue (dark-grey) time 
region starts T = 6 ms after the begining of the pre- 
vious time region, and presents a maximum amplitude 
of 1.3 N. We calculate PSDs for both time regions (fre- 
quency resolution Sf = 195 Hz) as shown in Fig. BJd). 
The PSDs shown in Fig. 0|d) correspond to the time re- 
gions of the same color shown in Fig. |4|c). Here, the 
PSDs are normalized by dividing the PSD by the peak 
PSD amplitude of the identified defect mode. It is evi- 
dent that the peak in the PSD spectrum corresponding 
to the time region with larger force amplitude presents 
a lower characteristic frequency (i.e., it is further into 
the gap) with respect to the peak representing the time 
region with lower force amplitudes. This is in agree- 
ment with the shift predicted by the parametric contin- 
uation analysis shown in Fig. HJa). The peak frequency 
of the PSD of the high force amplitude time region is 
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fdef — 13.28 kHz, and the peak frequency of the PSD 
of the low amplitude time region is fdef = 13.48 kHz, 
where fdef = 13.48 kHz is closer to the mean experimen- 
tally determined linear defect mode frequency (shown by 
the dashed line in Fig.QJd)). 
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FIG. 4: [Color online] (a) Numerical frequency continuation 
of the nonlinear defect modes corresponding to the experi- 
mental setup in Fig. QJa). (b) Numerically calculated spa- 
tial profile of the nonlinear localized mode with frequency 
fdef = 13.28 kHz. (c) Measured force-time history of a sen- 
sor at site n = 3, where a high amplitude, short width, force 
pulse is applied to the granular crystal, (d) Normalized PSD 
for the measured time regions of the same color in (c); closed 
and open diamonds correspond to the high and low amplitude 
time regions respectively. The vertical dashed line is the mean 
experimentally determined linear defect mode frequency. 



CONCLUSIONS 

We studied the response of statically compressed gran- 
ular crystals containing light-mass defects, and character- 
ized their near-linear spectra by applying continuous ex- 
citation. We demonstrated that such chains support lo- 
calized modes with frequencies above that of their acous- 
tic band cutoff, using approximate few-bead analytical 
calculations, numerics, and experiments. The number of 
supported localized modes depends on the number of de- 
fects, while their frequencies depend on the inter-defect 
distance, on the ratio jj of defect to regular masses (and 
the geometric/elastic properties of the beads), and on 
the static load. We also briefly described the nonlinear 
generalizations of such modes, departing from the near- 
linear regime, and showed a downshift of the correspond- 
ing defect mode frequencies with increasing amplitude. 
This study is important for understanding the interplay 
of disorder and nonlinearity in discrete systems, and the 
results reported may be relevant in the design of appli- 
cations involving vibrational energy trapping. 
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